#personalizar nuestra seccion
options(digits=3,papersize = "letter")
# https://stat.ethz.ch/R-manual/R-devel/library/base/html/options.html
# Opciones globales en graficas
par(mar=c(5.1,5,4.1,2.1),font=3,family="sans",bg="yellow")
Una empresa quiere introducir un programa de estímulos para su personal. Debe hacerlo lo más justo posible para evitar disgustos y un ambiente pesado entre los empleados. La información que ha recabado es de una muestra aleatoria de 200 empleados.
Antes de realizar el análisis estadístico más adecuado es conveniente revisar:
setwd("C:/Users/Emman/Documents/CCH-N} FES-A/fessaa/FES/7mo semestre/Temas Selectos de Estadistica")
getwd()
## [1] "C:/Users/Emman/Documents/CCH-N} FES-A/fessaa/FES/7mo semestre/Temas Selectos de Estadistica"
datos1 <- read.table("./DATA/Empleados.txt",header=TRUE,row.names=1)
head(datos1) # ver los primeros 6 renglones del conjunto de datos
## antiguedad horasextra sexo cursos incapacidad aptitudes escolaridad salario
## 1 11 125 1 4 9 121.9 2 23065
## 2 24 225 2 2 2 114.2 1 27180
## 3 17 115 2 3 5 134.1 1 34875
## 4 9 117 1 1 1 114.0 1 23685
## 5 15 26 1 2 0 151.4 2 33550
## 6 6 43 1 4 3 96.7 1 22635
## edad
## 1 44
## 2 50
## 3 48
## 4 53
## 5 62
## 6 45
tail(datos1) # ver los ultimos 6 renglones del conjunto de datos
## antiguedad horasextra sexo cursos incapacidad aptitudes escolaridad salario
## 195 5 173 2 2 8 122.1 2 29000
## 196 1 11 2 1 4 93.9 3 30200
## 197 7 129 2 3 6 104.6 1 34450
## 198 2 162 1 3 3 107.8 1 23550
## 199 2 5 2 1 1 101.7 1 22000
## 200 5 74 2 3 6 111.0 1 33000
## edad
## 195 29
## 196 30
## 197 35
## 198 25
## 199 26
## 200 34
dim(datos1)
## [1] 200 9
names(datos1)
## [1] "antiguedad" "horasextra" "sexo" "cursos" "incapacidad"
## [6] "aptitudes" "escolaridad" "salario" "edad"
str(datos1)
## 'data.frame': 200 obs. of 9 variables:
## $ antiguedad : int 11 24 17 9 15 6 4 2 17 17 ...
## $ horasextra : int 125 225 115 117 26 43 124 71 166 158 ...
## $ sexo : int 1 2 2 1 1 1 2 2 2 1 ...
## $ cursos : int 4 2 3 1 2 4 2 1 2 3 ...
## $ incapacidad: int 9 2 5 1 0 3 4 1 5 2 ...
## $ aptitudes : num 122 114 134 114 151 ...
## $ escolaridad: int 2 1 1 1 2 1 2 1 1 1 ...
## $ salario : int 23065 27180 34875 23685 33550 22635 19575 20430 18955 25595 ...
## $ edad : int 44 50 48 53 62 45 26 28 33 40 ...
mode(datos1)
## [1] "list"
# Exploremos algunas variables
datos1$antiguedad
## [1] 11 24 17 9 15 6 4 2 17 17 15 21 4 12 23 20 19 12 5 11 11 8 20 1 6
## [26] 18 21 7 21 27 20 11 11 3 16 2 12 16 9 15 3 17 17 23 6 1 7 4 22 18
## [51] 22 25 15 24 15 2 17 7 8 11 10 10 8 15 8 23 24 22 6 3 16 4 12 23 13
## [76] 11 4 3 9 12 3 16 23 2 18 3 7 25 2 17 22 6 27 14 12 24 14 14 11 4
## [101] 3 12 12 19 12 4 5 23 9 24 21 19 14 3 3 8 15 18 7 11 5 18 12 2 26
## [126] 26 11 11 0 7 19 5 26 1 8 3 3 14 11 22 12 17 20 11 14 22 2 14 5 19
## [151] 22 26 8 16 25 7 23 7 16 2 22 13 19 7 4 24 11 8 9 22 25 14 18 8 22
## [176] 8 13 27 27 3 2 23 9 15 17 5 26 27 12 1 3 0 16 3 5 1 7 2 2 5
datos1[,3] # da la columna 3
## [1] 1 2 2 1 1 1 2 2 2 1 2 2 1 1 1 1 2 2 2 2 1 2 1 1 2 1 1 1 2 2 2 1 1 1 1 2 1
## [38] 2 2 2 2 1 1 2 2 2 1 1 2 2 1 2 1 1 1 2 2 2 1 2 2 1 1 1 1 2 2 1 1 1 2 2 2 1
## [75] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 2 2 2 1 1 1 2 1 2 2 2 2 2 2 1 2 1 2 2 1 2 1
## [112] 2 1 1 2 2 1 1 2 2 1 1 1 1 1 1 2 2 2 1 1 2 2 2 2 2 1 1 1 2 1 2 2 1 1 2 1 2
## [149] 2 1 2 1 1 2 2 1 2 2 1 1 2 2 2 2 1 1 2 2 2 1 2 2 2 2 1 2 1 2 2 2 2 1 1 2 2
## [186] 2 2 1 2 2 2 2 2 1 2 2 2 1 2 2
datos1[3,] # da el renglon 3
## antiguedad horasextra sexo cursos incapacidad aptitudes escolaridad salario
## 3 17 115 2 3 5 134 1 34875
## edad
## 3 48
datos1$sexo # da la columna 3
## [1] 1 2 2 1 1 1 2 2 2 1 2 2 1 1 1 1 2 2 2 2 1 2 1 1 2 1 1 1 2 2 2 1 1 1 1 2 1
## [38] 2 2 2 2 1 1 2 2 2 1 1 2 2 1 2 1 1 1 2 2 2 1 2 2 1 1 1 1 2 2 1 1 1 2 2 2 1
## [75] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 2 2 2 1 1 1 2 1 2 2 2 2 2 2 1 2 1 2 2 1 2 1
## [112] 2 1 1 2 2 1 1 2 2 1 1 1 1 1 1 2 2 2 1 1 2 2 2 2 2 1 1 1 2 1 2 2 1 1 2 1 2
## [149] 2 1 2 1 1 2 2 1 2 2 1 1 2 2 2 2 1 1 2 2 2 1 2 2 2 2 1 2 1 2 2 2 2 1 1 2 2
## [186] 2 2 1 2 2 2 2 2 1 2 2 2 1 2 2
datos1[1:10,3:6] # para seleccionar una parte de los datos
## sexo cursos incapacidad aptitudes
## 1 1 4 9 121.9
## 2 2 2 2 114.2
## 3 2 3 5 134.1
## 4 1 1 1 114.0
## 5 1 2 0 151.4
## 6 1 4 3 96.7
## 7 2 2 4 98.4
## 8 2 1 1 110.1
## 9 2 2 5 102.0
## 10 1 3 2 101.0
# Si llamamos a la variable por su nombre...
#antiguedad # habiamos llamado datos1$antiguedad
#antiguedad*12
# MARCA ERROR AL TEJER!
# Solucion.....
attach(datos1) # con el attach ya puedo accesar a los nombres de las variables
antiguedad
## [1] 11 24 17 9 15 6 4 2 17 17 15 21 4 12 23 20 19 12 5 11 11 8 20 1 6
## [26] 18 21 7 21 27 20 11 11 3 16 2 12 16 9 15 3 17 17 23 6 1 7 4 22 18
## [51] 22 25 15 24 15 2 17 7 8 11 10 10 8 15 8 23 24 22 6 3 16 4 12 23 13
## [76] 11 4 3 9 12 3 16 23 2 18 3 7 25 2 17 22 6 27 14 12 24 14 14 11 4
## [101] 3 12 12 19 12 4 5 23 9 24 21 19 14 3 3 8 15 18 7 11 5 18 12 2 26
## [126] 26 11 11 0 7 19 5 26 1 8 3 3 14 11 22 12 17 20 11 14 22 2 14 5 19
## [151] 22 26 8 16 25 7 23 7 16 2 22 13 19 7 4 24 11 8 9 22 25 14 18 8 22
## [176] 8 13 27 27 3 2 23 9 15 17 5 26 27 12 1 3 0 16 3 5 1 7 2 2 5
antiguedad*12
## [1] 132 288 204 108 180 72 48 24 204 204 180 252 48 144 276 240 228 144
## [19] 60 132 132 96 240 12 72 216 252 84 252 324 240 132 132 36 192 24
## [37] 144 192 108 180 36 204 204 276 72 12 84 48 264 216 264 300 180 288
## [55] 180 24 204 84 96 132 120 120 96 180 96 276 288 264 72 36 192 48
## [73] 144 276 156 132 48 36 108 144 36 192 276 24 216 36 84 300 24 204
## [91] 264 72 324 168 144 288 168 168 132 48 36 144 144 228 144 48 60 276
## [109] 108 288 252 228 168 36 36 96 180 216 84 132 60 216 144 24 312 312
## [127] 132 132 0 84 228 60 312 12 96 36 36 168 132 264 144 204 240 132
## [145] 168 264 24 168 60 228 264 312 96 192 300 84 276 84 192 24 264 156
## [163] 228 84 48 288 132 96 108 264 300 168 216 96 264 96 156 324 324 36
## [181] 24 276 108 180 204 60 312 324 144 12 36 0 192 36 60 12 84 24
## [199] 24 60
# OJO: el problema con attach es que si manejo 2 o mas bases de datos distintas,
# y estas comparten el mismo nombre de las columnas, R se hace bolas y no
# sabra exactamente a que base de datos se refiere.
#detach(datos1) # para decirle a R que ya NO usare esa base de datos
# Datos faltantes?
library(naniar)
pct_miss(datos1)
## [1] 0
Antes de realizar el análisis estadístico más adecuado es conveniente revisar:
| Variable | Tipo | Escala de medición |
|---|---|---|
| Antigüedad | cuantitativa/discreta | razón |
| Hora extra | cuantitativa/discreta | razon |
| Sexo | cualitativa | nominal |
| Cursos | cuantitativa/discreta | razón |
| Incapacidad | cuantitativa/discreta | razon |
| Aptitudes | cuantitativa/continua | intervalo |
| Escolaridad | cualitativa | ordinal |
| Salarios | cuantitativa/continua | razón |
| Edad | cuantitativa/discreta | razón |
Basándonos en la distribución de frecuencias,
# Tabla de frecuencias
(tabla_sexo <- table(sexo)) # tabla de frecuencias absolutas
## sexo
## 1 2
## 81 119
(tablar_sexo <- prop.table(tabla_sexo)) # tabla de frecuencias relativas
## sexo
## 1 2
## 0.405 0.595
# otra forma de obtener las marginales:
(f.relativa_sexo<- c(sum(sexo==1)/length(sexo),sum(sexo==2)/length(sexo)))
## [1] 0.405 0.595
En resumen:
# Tabla de frecuencias completa
tabla.frec_sex<- matrix(cbind(tabla_sexo[1],tabla_sexo[2],
sum(tabla_sexo), tablar_sexo[1],tablar_sexo[2],sum(tablar_sexo)),
byrow=T,nrow = 2,ncol=3)
colnames(tabla.frec_sex)<- c("Mujer","Hombre","Totales")
rownames(tabla.frec_sex)<- c("fi","pi")
tabla.frec_sex
## Mujer Hombre Totales
## fi 81.000 119.000 200
## pi 0.405 0.595 1
Observamos que en esta empresa hay más empleados que empleadas. De los primeros hay 119 personas que constituyen el 59.5% del total.
Hacemos una gráfica para ver la distribución de los empleados (por género) en esta empresa:
# Diagrama de barras
barplot(tabla_sexo,names.arg=c("Mujer","Hombre"),col=c(3,32),
ylim=c(0,140))
abline(h=tabla_sexo,col=c("magenta","blue"),lwd=2) # dibuja las lineas horizontales
title("Diagrama de barras - Sexo \n (Frecuencias Absolutas)")
barplot(tablar_sexo,names.arg=c("M","H"),col=3:4,ylim=c(0,0.6))
title("Diagrama de barras - Sexo \n (Frecuencias Relativas)")
# Diagrama Circular
par(mfrow=c(1,2))
pie(tabla_sexo,labels=paste(c("M","H"),tabla_sexo,sep=": "),col=5:6,radius=1)
title("Diagrama de pastel - Sexo \n (Frecuencias absolutas)") # \n para escribir en otro renglon
pie(tablar_sexo,labels=paste(c("M","H"),tablar_sexo,sep=": "),col=15:16,radius=1)
title("Diagrama de pastel - Sexo \n (Frecuencias relativas)")
De la gráfica de barras y circular se observa que, en efecto, hay diferencias entre el número de empleados y empleadas que laboran en esta empresa. Vale la pena analizar estadísticamente si esto es así.
# Comparacion entre la proporcion de mujeres y la proporcion de hombres.
n_muj<- tabla_sexo[1]
n_hom<- tabla_sexo[2]
prop.test(x=c(n_muj, n_hom), n=c(200,200), conf.level=0.90)
##
## 2-sample test for equality of proportions with continuity correction
##
## data: c(n_muj, n_hom) out of c(200, 200)
## X-squared = 14, df = 1, p-value = 2e-04
## alternative hypothesis: two.sided
## 90 percent confidence interval:
## -0.276 -0.104
## sample estimates:
## prop 1 prop 2
## 0.405 0.595
De acuerdo con el intervalo de confianza para la diferencia de proporciones se concluye que SÍ hay diferencia entre el número de empleados y empleadas, de hecho \(p_{muj}-p_{hom}<0\); así que la proporción de mujeres si es estadísticamente inferior a la correspondiente de los hombres (como lo muestran las gráficas).
***NOTA: # ¿Tiene sentido la siguiente grafica?
hist(sexo,breaks=0:2) # eje vertical es la frecuencia absoluta
hist(sexo,breaks=0:2,probability=TRUE) # eje vertical es la frecuencia relativa
Sólo puede obtenerse la moda.
# install.packages(modeest)
library(modeest) # "genefilter" es una dependencia que NO esta disponible para R 3.6.0
# Otra paqueteria
# install.packages(modes)
# library(modes) # no disponible en R 4.0.0
# Variable nominal
mfv(sexo) # most frequent value
## [1] 2
mode(sexo) # NO!
## [1] "numeric"
# modes(sex) # library(modes)
summary(sexo) # no tiene interpretacion
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.00 1.00 2.00 1.59 2.00 2.00
Basándonos en la distribución de frecuencias,
# Tabla de frecuencias
(tabla_esc <- table(escolaridad)) # tabla de frecuencias absolutas
## escolaridad
## 0 1 2 3
## 34 117 47 2
(tablar_esc <- prop.table(tabla_esc)) # tabla de frecuencias relativas
## escolaridad
## 0 1 2 3
## 0.170 0.585 0.235 0.010
En resumen:
# Tabla de frecuencias completa
tabla.frec_esc<- matrix(cbind(tabla_esc[1],tabla_esc[2],
tabla_esc[3], tabla_esc[4],
sum(tabla_esc),
tablar_esc[1],tablar_esc[2],
tablar_esc[3],
tablar_esc[4], sum(tablar_esc)),
byrow=T,nrow = 2,ncol=5)
colnames(tabla.frec_esc)<- c("B","L(st)","L(ct)","P","Totales")
rownames(tabla.frec_esc)<- c("fi","pi")
tabla.frec_esc
## B L(st) L(ct) P Totales
## fi 34.00 117.000 47.000 2.00 200
## pi 0.17 0.585 0.235 0.01 1
Observamos que en esta empresa hay 117 emplead@s con licenciatura SIN título y representan el 58.5% del total. La diferencia entre los que sólo tienen el bachillerato y licenciatura CON título es pequeña (y veremos sí es significativa). Definitivamente casi no hay emplead@s con posgrado.
Hacemos una gráfica para ver la distribución de los empleados (por escolaridad) en esta empresa:
# Diagrama de barras
barplot(tabla_esc,names.arg=c("B","L(st)","L(ct)","P"),
col=c(3,32,2),ylim=c(0,140))
# Diagrama Circular
pie(tablar_esc,labels=paste(c("B","L(st)","L(ct)","P"),
tablar_esc,
sep=": "),col=rainbow(4),radius=1)
title("Diagrama de pastel - Escolaridad \n (Frecuencias relativas)")
De la gráfica de barras y circular se observa que, en efecto, hay diferencias entre el nivel educativo de los empleado(a)s. Vale la pena analizar estadísticamente si esto es así.
# Comparacion entre la proporcion de mujeres y la proporcion de hombres.
n_bach<- tabla_esc[1]
n_lct<- tabla_esc[3]
prop.test(x=c(n_bach, n_lct), n=c(200,200), conf.level=0.90)
##
## 2-sample test for equality of proportions with continuity correction
##
## data: c(n_bach, n_lct) out of c(200, 200)
## X-squared = 2, df = 1, p-value = 0.1
## alternative hypothesis: two.sided
## 90 percent confidence interval:
## -0.13588 0.00588
## sample estimates:
## prop 1 prop 2
## 0.170 0.235
Basándonos en el valor-p, NO se rechaza H0 por lo que \(p_{bach}-p_{lct}=0\) y por consiguiente no hay evidencia estadística suficiente para decir que la proporción de emplead@s con licenciatura CON título es diferente a la correspondiente con bachillerato.
Sólo tienen sentido: moda y mediana (o cualquier porcentil)
mfv(escolaridad)
## [1] 1
Los emplead@s con licenciatura SIN titulo son los que más laboran en esta empresa
median(escolaridad)
## [1] 1
quantile(escolaridad,prob=c(0.10,0.25,0.5,0.75))
## 10% 25% 50% 75%
## 0 1 1 1
El 10% de los emplead@s tienen a lo más bachillerato. Así que el 90% restante tienen educación superior. WOW!
El 75% de los emplead@s tienen a lo más licenciatura SIN título, por lo que el 25% restante ya tienen el título o poseen un posgrado.
Se trata de una variable discreta. Así que podemos presentar su resumen mediante una distribución de frecuencias como la correspondiente a una variable cualitativa.
# Tabla de frecuencias
(tabla_cursos <- table(cursos)) # tabla de frecuencias absolutas
## cursos
## 0 1 2 3 4 5 6 9
## 27 41 43 41 27 19 1 1
(tablar_cursos <- prop.table(tabla_cursos)) # tabla de frecuencias relativas
## cursos
## 0 1 2 3 4 5 6 9
## 0.135 0.205 0.215 0.205 0.135 0.095 0.005 0.005
En resumen
# Tabla de frecuencias completa
tabla.frec_cursos<- matrix(cbind(tabla_cursos[1],tabla_cursos[2],
tabla_cursos[3], tabla_cursos[4],
tabla_cursos[5], tabla_cursos[6],
tabla_cursos[7], tabla_cursos[8],
sum(tabla_cursos),
tablar_cursos[1], tablar_cursos[2],
tablar_cursos[3], tablar_cursos[4],
tablar_cursos[5], tablar_cursos[6],
tablar_cursos[7], tablar_cursos[8],
sum(tablar_cursos)),
byrow=T,nrow = 2,ncol=9)
colnames(tabla.frec_cursos)<- c("0","1","2","3","4","5","6","9","Totales")
rownames(tabla.frec_cursos)<- c("fi","pi")
tabla.frec_cursos
## 0 1 2 3 4 5 6 9 Totales
## fi 27.000 41.000 43.000 41.000 27.000 19.000 1.000 1.000 200
## pi 0.135 0.205 0.215 0.205 0.135 0.095 0.005 0.005 1
Veamos la distribución del número de cursos para los emplead@s de esta empresa.
par(mfrow=c(1,2))
barplot(tabla_cursos) # frecuencias absolutas
title("Diagrama de barras (Cursos) \n (Frecuencias absolutas)")
barplot(tablar_cursos,col=c(10:20)) # frecuencias relativas
title("Diagrama de barras (Points) \n (Frecuencias relativas)")
¿Qúe pasa si a la ingenua graficamos un histograma?
hist(cursos,breaks=-1:9,plot=FALSE) # las clases son: [-1,0] (0,1],...,(8,9])
## $breaks
## [1] -1 0 1 2 3 4 5 6 7 8 9
##
## $counts
## [1] 27 41 43 41 27 19 1 0 0 1
##
## $density
## [1] 0.135 0.205 0.215 0.205 0.135 0.095 0.005 0.000 0.000 0.005
##
## $mids
## [1] -0.5 0.5 1.5 2.5 3.5 4.5 5.5 6.5 7.5 8.5
##
## $xname
## [1] "cursos"
##
## $equidist
## [1] TRUE
##
## attr(,"class")
## [1] "histogram"
# OJO: las clases (breaks) los defini asi para que se pareciera al diagrama
# de barras y poderlo comparar.
hist(cursos,breaks=-1:9,probability=TRUE)
El perfil en ambas gráficas es “idéntico”, sólo que en el segundo caso las barras están pegadas. ¿Será correcto?
Podemos pensar que los valores de esta variable discreta son como las categorías y, por consiguiente, el diagrama de pay sería adecuado:
par(mfrow=c(1,2))
pie(tabla_cursos,labels=paste(c("0","1","2","3","4","5","6","9"),
tabla_cursos,sep=": "),col=palette("Alphabet"),radius=1)
title("Diagrama de pay - Cursos \n (Frecuencias absolutas)")
pie(tablar_cursos,labels=paste(c("0","1","2","3","4","5","6","9"),
tablar_cursos,sep=": "),col=topo.colors(8),radius=1)
title("Diagrama de pay - Cursos \n (Frecuencias relativas)")
Por tratarse de una variable cuantitativa es posible representar sus valores mediante un diagrama de puntos (dotplot o strip chart).
dotchart(cursos,labels=row.names(datos1),cex=0.8,groups = as.factor(sexo),#forzar a sexo a ser categorias, si son, pero no lo sabe R
xlab="Número de cursos",main="Diagrama de puntos") # Cleveland dot plot
Mejor tratemos con un subconjunto:
datos2<- subset(datos1,cursos>4); datos2
## antiguedad horasextra sexo cursos incapacidad aptitudes escolaridad salario
## 32 11 140 1 5 7 94.8 0 27500
## 41 3 6 2 9 3 107.9 2 19000
## 46 1 20 2 5 4 104.7 2 27050
## 49 22 77 2 5 4 80.5 0 40900
## 66 23 74 2 5 4 110.2 1 40370
## 67 24 138 2 5 4 103.6 1 42450
## 73 12 0 2 6 8 104.8 1 35850
## 76 11 90 2 5 2 101.7 2 40700
## 80 12 136 2 5 4 113.2 2 41800
## 87 7 132 2 5 7 126.3 1 28215
## 118 18 44 1 5 5 96.2 1 36200
## 122 18 182 1 5 8 104.4 1 36100
## 131 19 13 1 5 3 119.1 2 38300
## 143 20 39 2 5 2 123.4 1 43945
## 147 2 90 1 5 6 103.5 1 23325
## 153 8 169 1 5 4 87.8 2 31625
## 166 24 37 1 5 4 99.1 1 47125
## 169 9 198 2 5 3 107.8 1 31570
## 177 13 129 1 5 4 104.2 1 38300
## 180 3 209 2 5 3 92.4 0 23600
## 190 1 70 2 5 4 96.4 0 20000
## edad
## 32 32
## 41 23
## 46 31
## 49 63
## 66 65
## 67 63
## 73 61
## 76 42
## 80 42
## 87 31
## 118 51
## 122 46
## 131 40
## 143 62
## 147 25
## 153 34
## 166 64
## 169 44
## 177 38
## 180 22
## 190 20
dotchart(datos2$cursos,labels=row.names(datos2),cex=0.8,xlab="Número de cursos",main="Diagrama de puntos") # Cleveland dot plot
¿Y si mejor seleccionamos al azar?
#si queremos que sea lo mismo poner una semilla
datos3<- sample(datos1$cursos,size=60,replace = F); datos3 #replace F es igual a sin remplazo
## [1] 3 9 3 1 1 2 1 1 1 1 3 3 4 1 4 0 0 5 2 4 1 1 2 2 3 3 3 1 2 0 3 0 1 4 6 1 4 3
## [39] 1 3 5 2 4 2 5 2 2 1 2 1 1 1 5 2 2 3 3 4 4 5
dotchart(datos3,cex=0.8,xlab="Numero de cursos",main="Diagrama de puntos") # Cleveland dot plot
# Muestreo aleatorio de TODO el conjunto de datos
#sample(datos1,10,replace = F) # no es correcto!
datos4<- datos1[sample(nrow(datos1),30,replace = F),]; datos4
## antiguedad horasextra sexo cursos incapacidad aptitudes escolaridad salario
## 164 7 123 2 4 2 111.0 1 31575
## 56 2 169 2 1 9 118.4 1 15625
## 20 11 27 2 2 4 101.7 1 23960
## 31 20 99 2 0 4 96.9 2 37100
## 129 0 81 2 3 7 104.2 1 33900
## 30 27 157 2 2 5 113.2 2 39975
## 148 14 37 2 4 3 106.4 2 31325
## 106 4 200 1 2 3 100.9 2 32350
## 39 9 124 2 1 5 87.7 2 27200
## 168 8 100 2 2 6 107.8 2 40000
## 179 27 118 2 0 8 119.1 2 44000
## 149 5 173 2 2 4 110.2 1 34325
## 84 2 187 2 3 4 121.5 2 25150
## 122 18 182 1 5 8 104.4 1 36100
## 200 5 74 2 3 6 111.0 1 33000
## 69 6 68 1 2 5 87.6 1 39500
## 170 22 198 1 1 9 122.7 0 41750
## 189 12 209 2 2 3 117.5 1 31450
## 196 1 11 2 1 4 93.9 3 30200
## 159 16 105 1 2 5 89.7 1 39400
## 6 6 43 1 4 3 96.7 1 22635
## 90 17 118 2 1 3 98.1 0 40425
## 54 24 215 1 0 7 95.3 0 26180
## 124 2 217 1 2 0 123.1 1 34550
## 153 8 169 1 5 4 87.8 2 31625
## 107 5 182 2 1 7 124.1 1 34900
## 142 17 28 2 0 3 105.1 2 37675
## 94 14 193 1 3 4 122.0 1 42450
## 150 19 6 1 2 4 105.0 1 40500
## 123 12 48 1 4 3 111.4 1 28150
## edad
## 164 34
## 56 26
## 20 30
## 31 57
## 129 31
## 30 67
## 148 32
## 106 31
## 39 27
## 168 41
## 179 66
## 149 35
## 84 29
## 122 46
## 200 34
## 69 57
## 170 40
## 189 34
## 196 30
## 159 41
## 6 45
## 90 45
## 54 41
## 124 31
## 153 34
## 107 35
## 142 42
## 94 57
## 150 50
## 123 35
head(datos4)
## antiguedad horasextra sexo cursos incapacidad aptitudes escolaridad salario
## 164 7 123 2 4 2 111.0 1 31575
## 56 2 169 2 1 9 118.4 1 15625
## 20 11 27 2 2 4 101.7 1 23960
## 31 20 99 2 0 4 96.9 2 37100
## 129 0 81 2 3 7 104.2 1 33900
## 30 27 157 2 2 5 113.2 2 39975
## edad
## 164 34
## 56 26
## 20 30
## 31 57
## 129 31
## 30 67
dotchart(datos4$cursos,labels=row.names(datos4),cex=0.8,xlab="Numero de cursos",main="Diagrama de puntos") # Cleveland dot plot
Otro tipo de gráfica de puntos conocida como diagrama de puntos
#install.packages("BHH2")
library(BHH2)
dotPlot(na.omit(datos4$cursos),main="Diagrama de puntos (muestra aleatoria)",xlab="Número de puntos",pch=16)
Analicemos el diagrama de tallo y hojas. Recordemos que en este diagrama se divide a cada observación en dos partes: tallo y hoja.
El argumento “scale” permite modificar la longitud del diagrama y con ello tener este único decimal. Por ejemplo:
stem(datos1$cursos, scale=0.1)
##
## The decimal point is 1 digit(s) to the right of the |
##
## 0 | 00000000000000000000000000011111111111111111111111111111111111111111+99
## 0 | 555555555555555555569
Los valores que se leen son: 0.0,0.1,0.5,0.6,0.9
stem(datos1$cursos, scale=1) #no se usa mucho en ciendia de datos
##
## The decimal point is at the |
##
## 0 | 000000000000000000000000000
## 0 |
## 1 | 00000000000000000000000000000000000000000
## 1 |
## 2 | 0000000000000000000000000000000000000000000
## 2 |
## 3 | 00000000000000000000000000000000000000000
## 3 |
## 4 | 000000000000000000000000000
## 4 |
## 5 | 0000000000000000000
## 5 |
## 6 | 0
## 6 |
## 7 |
## 7 |
## 8 |
## 8 |
## 9 | 0
Los valores que se leen son: 0,1,2,3,4,5,6,9
COMENTARIO: En el primer caso, el valor exacto se pierde. Sin embargo, esto no es relevante para este tipo de gráfica/tabla. Sólo interesa identificar alrededor de cual valor se concentran más las observaciones y el perfil que muestra esta distribución.
Finalmente analizamos el diagrama de caja y brazos:
boxplot(cursos,horizontal=T,col="yellow")
Se observa que hay una persona con más de 8 cursos y es un outlier. Claramente la distribución de cursos es asimétrica hacia la derecha. El 50% de los emplead@s ha tomado a lo más 2 cursos; mientras que el 75% ha llevado hasta 3 cursos. Veamos de quién se trata para premiarlo.
datos1[datos1$cursos>8,]
## antiguedad horasextra sexo cursos incapacidad aptitudes escolaridad salario
## 41 3 6 2 9 3 108 2 19000
## edad
## 41 23
datos1[datos1$cursos>8,]$salario
## [1] 19000
library(fdth)
##
## Attaching package: 'fdth'
## The following object is masked from 'package:modeest':
##
## mfv
## The following objects are masked from 'package:stats':
##
## sd, var
mfv(cursos) #Usa libreria obligatoriamente
## [1] 2
mean(cursos)
## [1] 2.34
quantile(cursos, probs = c(0.25,0.50,0.75))
## 25% 50% 75%
## 1 2 3
IQR(cursos)
## [1] 2
var(cursos); sd(cursos)
## [1] 2.57
## [1] 1.6
mean(abs(na.omit(cursos)-mean(na.omit(cursos)))) # error medio (media)
## [1] 1.33
mean(abs(cursos-median(cursos))) # error medio (mediana)
## [1] 1.29
mean(abs(cursos-mfv(cursos))) # error medio (moda) OJO!!!
## [1] 1.29
c.v<- sd(cursos)/mean(cursos); c.v # tenemos poca variavilidad
## [1] 0.685
A veces es difícil identificar si la distribución de frecuencias es asimétrica o simétrica mediante un histograma o un diagrama de barras. Para ello existen medidas estadísticas que nos permiten saber si el perfil es simétrico o no; además si es un perfil “picudo” o “plano”. Estas medidas son: coeficiente de asimetría y coeficiente de curtosis.
n<- length(cursos) # longitud del vector
mu3<- sum((cursos-mean(cursos))^3)/(n-1)
s3<- sd(cursos)^3
j3<- mu3/s3; j3
## [1] 0.478
# Usando R
library(agricolae) #se puede con libreria momentos
##
## Attaching package: 'agricolae'
## The following object is masked from 'package:modeest':
##
## skewness
(j3_a<- skewness(cursos))
## [1] 0.483
(j3_b<- moments::skewness(cursos))
## [1] 0.479
mu4<- sum((cursos-mean(cursos))^4)/(n-1)
s4<- sd(cursos)^4
j4<- mu4/s4; j4
## [1] 3.24
# Usando R
library(agricolae) #numero 3 distribucion normal
(j4_a<- kurtosis(cursos)+3) # da el exceso de curtosis
## [1] 3.3
(j4_b<- moments::kurtosis(aptitudes))
## [1] 3.39
O bien, podemos resumir parte de esta información en:
fivenum(cursos) # minimum, lower-hinge, median, upper-hinge, maximum)
## [1] 0 1 2 3 9
summary(cursos)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 1.00 2.00 2.34 3.00 9.00
COMENTARIO SOBRE “HINGES”: es similar a los cuartiles. En conjuntos grandes de datos, la diferencia entre ellos es muy pequeña. Según la definición de Tukey los “hinges” son:
Hinges inferior: es la mediana de la mitad inferior de todos los valores a la izquierda de la mediana de los datos.
Hinges superior: es la mediana de la mitad superior de todos los valores a la derecha de la mediana de los datos
Ejemplo 1: n=20 observaciones La mediana está en la posición 10.5. La mitad inferior comprende las posiciones 1 hasta 10; mientras que la superior desde la 11 hasta 20.
El “lower hinge” es la mediana de la mitad inferior y sería la posición 5.5.
El “upper hinge” es la mediana de la mitad superior y sería la posición 5.5 (comenzando desde la posición 11), por lo que a final de cuentas será la posición 15.5.
Ejemplo 2: n=21 observaciones
La mediana está en la posición 11. La mitad inferior comprende las posiciones 1 hasta 11; mientras que la superior desde la 11 hasta 21.
El “lower hinge” es la mediana de la mitad inferior y sería la posición 6.
El “upper hinge” es la mediana de la mitad superior y sería la posición 6 (comenzando desde la posición 11), por lo que a final de cuentas será la posición 16.
Se trata de una variable continua, así que si hacemos lo siguiente:
table(salario)
## salario
## 15625 18955 19000 19575 19775 20000 20430 22000 22025 22400 22635 23065 23100
## 1 1 1 1 1 2 1 1 1 1 1 1 2
## 23150 23325 23550 23600 23650 23685 23960 24450 24525 24650 24750 24850 25100
## 1 1 1 3 1 1 1 1 1 2 1 1 2
## 25150 25545 25595 26180 26250 26275 26500 26800 27050 27175 27180 27200 27300
## 1 1 1 1 1 1 3 1 1 1 1 1 1
## 27350 27500 27600 27940 27950 28090 28150 28215 28650 28750 29000 29100 29500
## 1 3 1 1 1 1 2 1 1 1 2 1 1
## 29900 30000 30200 30440 30495 30500 30600 30630 31200 31325 31400 31450 31570
## 1 1 1 1 1 1 1 1 1 1 1 2 1
## 31575 31625 32000 32350 32670 32775 33000 33110 33550 33900 33965 33975 34250
## 1 1 3 1 1 1 3 1 1 1 1 1 1
## 34325 34450 34550 34850 34875 34900 34975 35000 35200 35300 35600 35850 35875
## 2 1 2 1 2 1 1 2 2 1 1 2 1
## 36100 36195 36200 36475 37000 37100 37200 37225 37300 37375 37600 37675 37750
## 2 1 1 1 1 1 1 1 1 1 2 1 2
## 38170 38300 38500 38790 38900 38950 39275 39375 39400 39500 39750 39975 40000
## 1 3 2 1 1 1 1 1 1 3 1 1 1
## 40050 40370 40400 40425 40500 40700 40900 41150 41575 41600 41650 41705 41750
## 1 1 1 1 3 1 1 1 1 2 1 1 1
## 41800 41900 41910 42000 42250 42450 42850 43175 43650 43675 43780 43800 43850
## 3 1 1 1 1 2 1 1 2 1 1 2 2
## 43945 44000 45000 45100 45375 45800 46100 46550 46900 47025 47100 47125 47200
## 1 2 1 1 2 1 1 1 1 1 1 1 1
## 47300 47950
## 1 1
barplot(salario)
COMENTARIO: Se observa que no es util para este tipo de variable representarlo asi, por eso conviene hacer una tabla de frecuencias usando INTERVALOS DE CLASE.
En tenemos dos opciones para construir una tabla de frecuencias para datos agrupados.
Sugerencia: revisar primero los valores mínimo y máximo
#imagen del logo en sea, guardado en el directorio de trabajo
range(salario) # minimo y maximo
## [1] 15625 47950
Con ellos, ahora sí definimos los breaks del histograma:
hist_sal <- hist(datos1$salario,breaks=seq(15600,49000,2500),plot=F)
summary(hist_sal)
## Length Class Mode
## breaks 14 -none- numeric
## counts 13 -none- numeric
## density 13 -none- numeric
## mids 13 -none- numeric
## xname 1 -none- character
## equidist 1 -none- logical
Construimos la correspondiente tabla de frecuencias
n <- length(hist_sal$breaks)
tab_sal <- cbind(hist_sal$breaks[-n],hist_sal$breaks[-1],
hist_sal$counts,
hist_sal$counts/sum(hist_sal$counts),
cumsum(hist_sal$counts),
cumsum(hist_sal$counts/sum(hist_sal$counts)))
dimnames(tab_sal)[[2]]<-c("Linf","Lsup","f","fr","F","Fr")
print(tab_sal)
## Linf Lsup f fr F Fr
## [1,] 15600 18100 1 0.005 1 0.005
## [2,] 18100 20600 7 0.035 8 0.040
## [3,] 20600 23100 7 0.035 15 0.075
## [4,] 23100 25600 20 0.100 35 0.175
## [5,] 25600 28100 20 0.100 55 0.275
## [6,] 28100 30600 16 0.080 71 0.355
## [7,] 30600 33100 18 0.090 89 0.445
## [8,] 33100 35600 22 0.110 111 0.555
## [9,] 35600 38100 19 0.095 130 0.650
## [10,] 38100 40600 25 0.125 155 0.775
## [11,] 40600 43100 19 0.095 174 0.870
## [12,] 43100 45600 16 0.080 190 0.950
## [13,] 45600 48100 10 0.050 200 1.000
library(fdth)
(tabla_cont<- fdt(salario,breaks="Scott",start=15600,end=49000,h=2500,right=T))
## Class limits f rf rf(%) cf cf(%)
## (15600,18100] 1 0.00 0.5 1 0.5
## (18100,20600] 7 0.04 3.5 8 4.0
## (20600,23100] 7 0.04 3.5 15 7.5
## (23100,25600] 20 0.10 10.0 35 17.5
## (25600,28100] 20 0.10 10.0 55 27.5
## (28100,30600] 16 0.08 8.0 71 35.5
## (30600,33100] 18 0.09 9.0 89 44.5
## (33100,35600] 22 0.11 11.0 111 55.5
## (35600,38100] 19 0.10 9.5 130 65.0
## (38100,40600] 25 0.12 12.5 155 77.5
## (40600,43100] 19 0.10 9.5 174 87.0
## (43100,45600] 16 0.08 8.0 190 95.0
## (45600,48100] 10 0.05 5.0 200 100.0
(tabla_cont<- fdt(salario,breaks="Sturges",right=T))
## Class limits f rf rf(%) cf cf(%)
## (15468.75,19131.06] 3 0.01 1.5 3 1.5
## (19131.06,22793.36] 9 0.04 4.5 12 6.0
## (22793.36,26455.67] 26 0.13 13.0 38 19.0
## (26455.67,30117.97] 28 0.14 14.0 66 33.0
## (30117.97,33780.28] 25 0.12 12.5 91 45.5
## (33780.28,37442.58] 34 0.17 17.0 125 62.5
## (37442.58,41104.89] 32 0.16 16.0 157 78.5
## (41104.89,44767.19] 29 0.14 14.5 186 93.0
## (44767.19,48429.5] 14 0.07 7.0 200 100.0
(tabla_cont<- fdt(salario,breaks="Sturges",start=15600,end=49000,h=2500,right=T))
## Class limits f rf rf(%) cf cf(%)
## (15600,18100] 1 0.00 0.5 1 0.5
## (18100,20600] 7 0.04 3.5 8 4.0
## (20600,23100] 7 0.04 3.5 15 7.5
## (23100,25600] 20 0.10 10.0 35 17.5
## (25600,28100] 20 0.10 10.0 55 27.5
## (28100,30600] 16 0.08 8.0 71 35.5
## (30600,33100] 18 0.09 9.0 89 44.5
## (33100,35600] 22 0.11 11.0 111 55.5
## (35600,38100] 19 0.10 9.5 130 65.0
## (38100,40600] 25 0.12 12.5 155 77.5
## (40600,43100] 19 0.10 9.5 174 87.0
## (43100,45600] 16 0.08 8.0 190 95.0
## (45600,48100] 10 0.05 5.0 200 100.0
(tabla_cont<- fdt(salario,breaks="FD",right=T))
## Class limits f rf rf(%) cf cf(%)
## (15468.75,19588.84] 4 0.02 2.0 4 2.0
## (19588.84,23708.94] 19 0.10 9.5 23 11.5
## (23708.94,27829.03] 29 0.14 14.5 52 26.0
## (27829.03,31949.12] 28 0.14 14.0 80 40.0
## (31949.12,36069.22] 34 0.17 17.0 114 57.0
## (36069.22,40189.31] 35 0.17 17.5 149 74.5
## (40189.31,44309.41] 37 0.18 18.5 186 93.0
## (44309.41,48429.5] 14 0.07 7.0 200 100.0
Hacemos el histograma:
par(mfrow=c(1,2))
hist(salario,col=rainbow(6), prob=T)
hist(salario,col=rainbow(6), prob=F)
El histograma muestra el perfil de la distribución de probabilidad que mejor describe a la variable. Para ello, resulta conveniente dibujar el polígono de frecuencias. De hecho,
hist(salario,probability=T,col="grey50") # histograma cuyo eje vertical esta en escala
# de "probabilidades".
lines(density(salario),lwd=2,col="blue") # curva de densidad de probabilidad
par(mfrow=c(1,2))
library(MASS)
frequency.polygon(salario,nclass = nclass.freq(salario))
## Warning in hist.default(x, nclass, probability = TRUE, plot = FALSE, ...):
## argument 'probability' is not made use of
library(UsingR)
## Loading required package: HistData
## Loading required package: Hmisc
## Loading required package: lattice
## Loading required package: survival
## Loading required package: Formula
## Loading required package: ggplot2
##
## Attaching package: 'Hmisc'
## The following objects are masked from 'package:base':
##
## format.pval, units
##
## Attaching package: 'UsingR'
## The following object is masked from 'package:survival':
##
## cancer
simple.freqpoly(salario,col=heat.colors(6))
NOTA: si se requieren intervalos desiguales para observar ciertos detalles de interés
hist_uneq<- hist(salario,breaks = quantile(salario, 0:10 / 10),col="gray") # usa los cuantiles
hist(salario,breaks=c(15600,seq(27000,32000,800),37690,49000)) # intervalos en puntos arbitrarios
Aprovechando la tabla de frecuencias con intervalos de clase y su histograma, vale la pena revisar la siguiente gráfica: ojiva
library(agricolae)
ogive.freq(hist_sal,main="Ojiva (salario)",xlab="Límite superior de
la clase",ylab="Pi (%)")
## Límite superior de\n la clase RCF
## 1 15600 0.000
## 2 18100 0.005
## 3 20600 0.040
## 4 23100 0.075
## 5 25600 0.175
## 6 28100 0.275
## 7 30600 0.355
## 8 33100 0.445
## 9 35600 0.555
## 10 38100 0.650
## 11 40600 0.775
## 12 43100 0.870
## 13 45600 0.950
## 14 48100 1.000
## 15 50600 1.000
# "ogive.freq" requiere al histograma como objeto
¿Cómo se interpreta esta gráfica?
# Haremos una grafica interactiva para Markdown.
# Primero debemos crear un data.frame con las columnas de la tabla de # frecuencias (requisito de la funcion "plot_ly")
library(plotly)
## Registered S3 method overwritten by 'httr':
## method from
## print.response rmutil
##
## Attaching package: 'plotly'
## The following object is masked from 'package:Hmisc':
##
## subplot
## The following object is masked from 'package:ggplot2':
##
## last_plot
## The following object is masked from 'package:MASS':
##
## select
## The following object is masked from 'package:stats':
##
## filter
## The following object is masked from 'package:graphics':
##
## layout
h<- data.frame(cbind(tab_sal[,2],tab_sal[,6])); h
## X1 X2
## 1 18100 0.005
## 2 20600 0.040
## 3 23100 0.075
## 4 25600 0.175
## 5 28100 0.275
## 6 30600 0.355
## 7 33100 0.445
## 8 35600 0.555
## 9 38100 0.650
## 10 40600 0.775
## 11 43100 0.870
## 12 45600 0.950
## 13 48100 1.000
plot_ly(data=h,x= ~tab_sal[,2],y= ~tab_sal[,6],
marker = list(size = 10, color="red"),
type = 'scatter', mode = 'lines')%>%
layout(title = 'Ojiva (Salario)',
xaxis = list(title = 'Límite superior'),
yaxis = list(title = 'Pi'),
shapes = list(
# linea vertical
list(type = "line", x0 = 0, x1 = 0,
y0 = 0, y1 = 1, yref = "paper"),
# linea horizontal
list(type = "line", x0 = 0, x1 = 50000,
y0 = 1, y1 = 1, xref = "paper")))
## A marker object has been specified, but markers is not in the mode
## Adding markers to the mode...
Veamos otras gráficas que ya habíamos realizado en el caso discreto.
dotPlot(salario,pch=16,col=5) # libreria(BHH2)
stem(salario)
##
## The decimal point is 3 digit(s) to the right of the |
##
## 14 | 6
## 16 |
## 18 | 0068
## 20 | 004
## 22 | 004611123666677
## 24 | 055778911256
## 26 | 233555812223455569
## 28 | 012227800159
## 30 | 024556623455666
## 32 | 000478000169
## 34 | 0033356699990002236999
## 36 | 1122501223466788
## 38 | 2333558903445558
## 40 | 00144455579266677888899
## 42 | 035592777888999
## 44 | 0001448
## 46 | 16901123
## 48 | 0
boxplot(salario,horizontal = T,col="lightgreen")
# Mejor...
library(plotly)
plot_ly(datos1,x = ~salario, type = "box")
¿Qué resultado muestra lo siguiente?
mfv(salario)
## [1] 27500
¿Tiene sentido para una variable continua?
Nótese que si las observaciones se redondean es posible que encontremos un valor que se repita varias veces. El programa R lo que hace es redondear hacia abajo (floor) y con ello localiza la moda. Esto no es adecuado.
Si conocemos la distribución de probabilidad para la variable “salario”, entonces basta con derivarla para tener el valor que maximiza su probabilidad. ¿Y si no, qué hacemos?
Respuesta: usamos su tabla de frecuencias
# Habiamos visto que
tabla_cont # tabla de frecuencias para datos agrupados por intervalos de clase
## Class limits f rf rf(%) cf cf(%)
## (15468.75,19588.84] 4 0.02 2.0 4 2.0
## (19588.84,23708.94] 19 0.10 9.5 23 11.5
## (23708.94,27829.03] 29 0.14 14.5 52 26.0
## (27829.03,31949.12] 28 0.14 14.0 80 40.0
## (31949.12,36069.22] 34 0.17 17.0 114 57.0
## (36069.22,40189.31] 35 0.17 17.5 149 74.5
## (40189.31,44309.41] 37 0.18 18.5 186 93.0
## (44309.41,48429.5] 14 0.07 7.0 200 100.0
str(tabla_cont) # para ver como debemos llamar a cada elemento
## List of 2
## $ table :'data.frame': 8 obs. of 6 variables:
## ..$ Class limits: Factor w/ 8 levels "(15468.75,19588.84]",..: 1 2 3 4 5 6 7 8
## ..$ f : int [1:8] 4 19 29 28 34 35 37 14
## ..$ rf : num [1:8] 0.02 0.095 0.145 0.14 0.17 0.175 0.185 0.07
## ..$ rf(%) : num [1:8] 2 9.5 14.5 14 17 17.5 18.5 7
## ..$ cf : num [1:8] 4 23 52 80 114 149 186 200
## ..$ cf(%) : num [1:8] 2 11.5 26 40 57 74.5 93 100
## $ breaks: Named num [1:4] 15469 48430 4120 1
## ..- attr(*, "names")= chr [1:4] "start" "end" "h" "right"
## - attr(*, "class")= chr [1:3] "fdt.default" "fdt" "list"
which.max(tabla_cont$table$f) # indica el renglon cuya frecuencia absoluta, f, es maxima
## [1] 7
tabla_cont$table$`Class limits`[which.max(tabla_cont$table$f)] # indica la CLASE MODAL
## [1] (40189.31,44309.41]
## 8 Levels: (15468.75,19588.84] (19588.84,23708.94] ... (44309.41,48429.5]
Comparamos los resultados y claramente son ¡MUY DISTINTOS!
mean(salario)
## [1] 34009
quantile(salario, probs = c(0.25,0.50,0.75))
## 25% 50% 75%
## 27500 34700 40378
IQR(salario)
## [1] 12878
var(salario); sd(salario)
## [1] 57075260
## [1] 7555
mean(abs(na.omit(salario)-mean(na.omit(salario)))) # error medio (media)
## [1] 6422
mean(abs(salario-median(salario))) # error medio (mediana)
## [1] 6398
mean(abs(salario-mfv(salario))) # error medio (moda) OJO!!!
## [1] 8245
c.v<- sd(salario)/mean(salario); c.v
## [1] 0.222
n<- length(salario) # longitud del vector
mu3<- sum((salario-mean(salario))^3)/(n-1)
s3<- sd(salario)^3
j3<- mu3/s3; j3
## [1] -0.143
# Usando R
library(agricolae)
(j3_a<- skewness(salario))
## [1] -0.144
(j3_b<- moments::skewness(salario))
## [1] -0.143
mu4<- sum((salario-mean(salario))^4)/(n-1)
s4<- sd(salario)^4
j4<- mu4/s4; j4
## [1] 2.03
# Usando R
library(agricolae)
(j4_a<- kurtosis(salario)+3) # da el exceso de curtosis
## [1] 2.05
(j4_b<- moments::kurtosis(salario))
## [1] 2.04
O bien, podemos resumir parte de esta información en:
fivenum(salario) # minimum, lower-hinge, median, upper-hinge, maximum)
## [1] 15625 27500 34700 40385 47950
summary(salario)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 15625 27500 34700 34009 40378 47950
Para comparar el comportamiento de estas variables conviene usar un diagrama de caja-brazos. Este diagrama se usa para la variable cuantitativa, la cualitativa solo sirve para separar el análisis por casos.
Escogemos las variables escolaridad y salario.
El diagrama de caja-brazos de salario se muestra en la sección 2.4.2.1. Es claro que hay asimetría hacia la izquierda (\(\gamma_{3}=-0.143\)) pero con esta gráfica aún no podemos saber si la escolaridad juega un papel relevante. Para ello hacemos el siguiente diagrama:
# Ahora, analicemos esta variable por escolaridad
sal_bach<- datos1[escolaridad==0,]$salario
boxplot(sal_bach,horizontal = T)
boxplot(salario~escolaridad,col="yellow")
# Otra manera
par(mfrow=c(1,1)) # ventana original
library(paletteer)
library(ggthemes)
boxplot(salario[escolaridad=="0"],salario[escolaridad=="1"],
salario[escolaridad=="2"],salario[escolaridad=="3"],
col=paletteer_c("ggthemes::Sunset-Sunrise Diverging", 4),names=c("bachillerato","licenciatura sin título","licenciatura con título","posgrado"))
# Otra manera
library(ggplot2)
ggplot(datos1,
aes(x = as.factor(escolaridad),
y = salario)) +
geom_boxplot() +
labs(title = "Salary distribution by rank")
El salario no es el mismo para todos los niveles de escolaridad. Las medianas de los niveles Lst, Lct y posgrado, se parecen. ¿Qué tal si hacemos una prueba de hipótesis para ver esto?
# Prueba NO PARAMETRICA para comparacion de medianas
wilcox.test(salario[escolaridad=="1"],
salario[escolaridad=="2"], alternative= "two.sided", mu=0, paired= F, conf.level= 0.97)
##
## Wilcoxon rank sum test with continuity correction
##
## data: salario[escolaridad == "1"] and salario[escolaridad == "2"]
## W = 2712, p-value = 0.9
## alternative hypothesis: true location shift is not equal to 0
De acuerdo con el valor-p se concluye que no hay evidencia estadística suficiente para pensar que la diferencia de medianas no es igual a cero, en otras palabras, la mediana de los salarios para los emplead@s con licenciatura sin título es igual a la correspondiente con licenciatura con título (¡ganas lo mismo con el “papel” que sin él!)
Y si comparamos a los que tienen posgrado:
# Prueba NO PARAMETRICA para comparacion de medianas
wilcox.test(salario[escolaridad=="1"],
salario[escolaridad=="3"], alternative= "two.sided", mu=0, paired= F, conf.level= 0.97)
##
## Wilcoxon rank sum test with continuity correction
##
## data: salario[escolaridad == "1"] and salario[escolaridad == "3"]
## W = 92, p-value = 0.6
## alternative hypothesis: true location shift is not equal to 0
¡Ufff! A la empresa no le importa que tengamos un posgrado porque me pagará igual que a la persona que sólo tiene el 100% de créditos de la carrera.
Veamos esto:
library(ggplot2)
# Distribucion del salario
Nivel.educacion<- as.factor(escolaridad)
num.empleado<- seq(1:200)
ggplot(datos1, aes(x = num.empleado,
y = salario,
color=Nivel.educacion)) +
geom_point() +
labs(title = "Salario de acuerdo a la escolaridad")
# O bien,
colors <- c("red", "green", "blue", "orange")
plot(num.empleado,salario,col = colors[Nivel.educacion],
pch = as.numeric(Nivel.educacion))
legend("bottomright", legend = c("B", "Lst","Lct","P"),
lwd = 3, col = c("red", "green", "blue", "orange"))
De acuerdo a las gráficas anteriores, no se observan grupos salariales.
Más sobre el salario:
plot(antiguedad, salario, col=colors[Nivel.educacion], pch = 16)
legend("bottomright", legend = c("B", "Lst","Lct","P"),
lwd = 3, col = c("red", "green", "blue", "orange"))
# calculate mean salary for each rank
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:Hmisc':
##
## src, summarize
## The following object is masked from 'package:MASS':
##
## select
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
plotdata <- datos1 %>%
group_by(escolaridad) %>%
summarize(sal_prom = mean(salario))
# plot mean salaries
ggplot(plotdata,
aes(x = escolaridad,
y = sal_prom)) +
geom_bar(stat = "identity")
library(scales)
ggplot(plotdata,
aes(x = factor(escolaridad,
labels = c("Bachillerato",
"Licenciatura sin título",
"Licenciatura con título","Posgrado")),
y = sal_prom)) +
geom_bar(stat = "identity",
fill = "cornflowerblue") +
geom_text(aes(label = dollar(sal_prom)),
vjust = -0.25) +
scale_y_continuous(breaks = seq(0, 50000, 10000),
label = dollar) +
labs(title = "Salario promedio por nivel educativo",
subtitle = "Empresa X",
x = "",
y = "")
Para estudiar la relación entre una variable que agrupa (cualitativa) y una cuantitativa podemos usar: strip plot
# Graficamos la distribucion de los salarios por nivel educativo usando strip plots
ggplot(datos1,
aes(y = escolaridad,
x = salario)) +
geom_point() +
labs(title = "Distribución de los salarios por nivel educativo")
Es difícil comprender lo que sucede, mejor:
library(scales)
ggplot(datos1,
aes(y = factor(escolaridad,labels=c("Bachillerato","Licenciatura\n sin título","Licenciatura\ncon título","Posgrado")),
x = salario,color = Nivel.educacion)) +
geom_jitter() +
scale_x_continuous(label = dollar) +
labs(title = "Distribución de los salarios por nivel educativo")+
labs(y= "Educación", x = "Salario (USD)")+
theme_minimal()
round(cor(datos1[c(1,2,4:6,8:9)]),digits = 2)
## antiguedad horasextra cursos incapacidad aptitudes salario edad
## antiguedad 1.00 0.07 0.00 0.05 0.07 0.63 0.76
## horasextra 0.07 1.00 -0.02 0.10 0.14 0.06 0.05
## cursos 0.00 -0.02 1.00 -0.13 -0.10 0.01 0.05
## incapacidad 0.05 0.10 -0.13 1.00 -0.03 0.00 0.04
## aptitudes 0.07 0.14 -0.10 -0.03 1.00 0.16 0.09
## salario 0.63 0.06 0.01 0.00 0.16 1.00 0.75
## edad 0.76 0.05 0.05 0.04 0.09 0.75 1.00
Observamos que las relaciones más débiles se dan entre: antiguedad-horas extra; antiguedad-cursos; antiguedad-incapacidad; antiguedad-aptitudes, horasextra-cursos, cursos-salario, cursos-edad, incapacidad-aptitudes, incapacidad-edad, edad-aptitudes.
De hecho, estas parejas de variables tienen coeficientes cercanos a cero. Así que podríamos sospechar que no hay relación lineal entre ellas. Para verificar esto, hagamos una prueba de hipótesis:
cor.test(antiguedad, horasextra, alternative="two.sided", method="pearson")
##
## Pearson's product-moment correlation
##
## data: antiguedad and horasextra
## t = 1, df = 198, p-value = 0.3
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.065 0.211
## sample estimates:
## cor
## 0.0744
cor.test(antiguedad, cursos, alternative="two.sided", method="pearson")
##
## Pearson's product-moment correlation
##
## data: antiguedad and cursos
## t = -0.04, df = 198, p-value = 1
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.141 0.136
## sample estimates:
## cor
## -0.00266
cor.test(antiguedad, incapacidad, alternative="two.sided", method="pearson")
##
## Pearson's product-moment correlation
##
## data: antiguedad and incapacidad
## t = 0.7, df = 198, p-value = 0.5
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.0895 0.1873
## sample estimates:
## cor
## 0.0499
cor.test(antiguedad, aptitudes, alternative="two.sided", method="pearson")
##
## Pearson's product-moment correlation
##
## data: antiguedad and aptitudes
## t = 1, df = 198, p-value = 0.3
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.0647 0.2112
## sample estimates:
## cor
## 0.0747
cor.test(cursos, horasextra, alternative="two.sided", method="pearson")
##
## Pearson's product-moment correlation
##
## data: cursos and horasextra
## t = -0.3, df = 198, p-value = 0.8
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.160 0.118
## sample estimates:
## cor
## -0.0213
Gráficamente…
plot(datos1[c(1,2,4:6,8:9)]) # matriz de graficas de dispersion
Comparamos tanto el salario como las aptitudes:
cv_aptitudes<- sd(aptitudes)/mean(aptitudes); cv_aptitudes
## [1] 0.118
cv_salario<- sd(salario)/mean(salario); cv_salario
## [1] 0.222
¿Dónde hay más variabilidad en cuanto al salario de acuerdo con la escolaridad?
bachillerato<- datos1[escolaridad=="0",]
lic.st<- datos1[escolaridad=="1",]
lic.ct<- datos1[escolaridad=="2",]
posgrado<- datos1[escolaridad=="3",]
(cv_bachillerato<- sd(bachillerato$salario)/mean(bachillerato$salario))
## [1] 0.235
(cv_lic.st<- sd(lic.st$salario)/mean(lic.st$salario))
## [1] 0.222